clear all;clc;close all;
w=0.2;OD=10;power=20;P0=1.0;
cs = CubeSolver('nGrid',100,'lengthD',1.0,'OD',OD,'power',power,'P0',P0);
beamshape=@(x,y,z) exp(-((x-0.5).*(x-0.5) + (y-0.5).*(y-0.5))/2/w^2).*boolean(x>0.5);
beamshape=@(x,y,z) exp(-((x-0.25).*(x-0.25) + (y-0.25).*(y-0.25))/2/w^2)+...
    exp(-((x-0.25).*(x-0.25) + (y-0.75).*(y-0.75))/2/w^2)+...
    exp(-((x-0.75).*(x-0.75) + (y-0.25).*(y-0.25))/2/w^2)+...
    exp(-((x-0.75).*(x-0.75) + (y-0.75).*(y-0.75))/2/w^2);
cs.dispPolarize(beamshape);